#!/usr/bin/perl  -w
use strict;
use Getopt::Long;

#the command line should be:  perl ripalign_analysis_new.pl 2010_4_26_17_57_49_ripalign.txt  -s sthe -c 1 for example
#and it wil generate a file like this TE_sthe_chrm1.tab 
#right now this only works for ndis,cglo,ntet,smac,sther,tter except for ncra



my @matrix;
my @temp;
#my $species;
#my $chrom;
my %rec;

#GetOptions("species=s" => \$species,
#	   "chrom=i" => \$chrom);
my $file = shift @ARGV;
open(DATA, $file) || die "cannot open file $file: $!\n";
my $outfile;

while(<DATA>)
{ 
	next if /^Repeat/; #leave the head of the file out	
	my @row=split(/\s+/,$_);
	if( ! defined $outfile ) {
	    # first time get to see a line with the TE family name
	    my @split = split(/_/,$row[0]);
	    $outfile = join("_",$split[-3], $split[-2], $split[-1]);
	}
	
	my $seqid=$row[8];
	my ($ca,$cc,$cg,$ct) = ($row[41],$row[48],$row[55],$row[62]);
	my @idinfor = split("_",$seqid);
	my ($sp,$chr,$st) = ($idinfor[0],$idinfor[2],$idinfor[3]);
	my $temp = $cc+$cg+$ct;

	my $ca_dom;	
	if($temp==0){$ca_dom=$ca;} #make sure that it won't be divided by zero
	else {
	    $ca_dom=$ca/($cc+$cg+$ct);  #calculate the CpA<->TpA dominance
	}
	# could also write it like this
	# my $ca_dom = ( $temp == 0 ) ? $ca : $ca/($cc+$cg+$ct);	
	push @{$rec{$sp}->{$chr}}, [$seqid,$sp,$chr,$st,$ca_dom];
}   

my %distance;


#this is to display the original matrix before sorting
#foreach my $row ( @rec) { 
 #       print join("\t",@$row),"\n";
#}

#printf"%-60s\t%-6s\t%3s\t%-20s\t%-10s\n","ID","SP","CHR","Start","CpA_Dom";
#this is to sort the matrix by the start position on the same chromosome

foreach my $species ( keys %rec ) {
    mkdir($species);
    my $output="$species/TE_".$outfile.".tab";
    open(OUT,">$output") || die $!;
    printf OUT "%-60s\t%-16s%12s\n","ID","Dist","CpA_Dom";
    
    foreach my $chrom ( keys %{ $rec{$species} } ) {

	my @rec = @{$rec{$species}->{$chrom}};
	if (@rec <= 1) {
	    next; # if there is only 1 entry for the chromosome, let's skip it
	}
	my (@id,@chrom,@start,@dom);

	foreach my $row (sort {$a->[3] <=> $b->[3]  } @rec) { 
#	printf "%-60s\t%-6s\t%3d\t%-20d\t%.2f\n",@$row[0],@$row[1],@$row[2],@$row[3],@$row[4];
	    push(@start,@$row[3]);
	    push(@id,@$row[0]);
	    push(@dom,@$row[4]);
	}

#printf"%-60s\t%16s%14s\n","ID","Dist","CpA_Dom";	
#this is to count the closest distance for each sequence at the chromosome
#	foreach my $each (@start)
#	{
#	print("$each\n");
#	}       
	for(my $j=0;$j<scalar(@start);$j++)
	{
	    my $dist;
	    my $post_dist; 
	    my $prev_dist; 
#for the first sequence,$dist=$post_dist
	    if($j==0) {
		$post_dist=$start[$j+1]-$start[$j];
		$dist = $post_dist;
		$distance{$id[$j]}=$dist;
	    } elsif($j==scalar(@start)-1) {
#for the last sequence,$dist=$prev_dist
		$prev_dist = $start[$j] - $start[$j-1];
		$dist      = $prev_dist;
		$distance{$id[$j]} = $dist;
	    } elsif($j>0 && $j<scalar(@start)-1){
		$post_dist=$start[$j+1]-$start[$j];	
		$prev_dist=$start[$j]-$start[$j-1];
		if ($post_dist<$prev_dist) {
		    $dist = $post_dist;
		    $distance{$id[$j]}=$dist;	
		} else {
		    $dist = $prev_dist;
		    $distance{$id[$j]}=$dist;
		}
	    }
	    printf OUT "%-60s\t%16s\t%6.2f\n", $id[$j],$dist,$dom[$j];
#   printf OUT "%-60s\t%16s\t%.2f\n", $id[$j],$dist,$dom[$j];
	}
    }
	close OUT;
}
#print("-------\n");
#foreach my $row ( @rec) { 
#       print join("\t",@$row[0],@$row[2],@$row[3]),"\n";
#}
#here we are going to make a new column which stores the shortest distance
